Analytical and Numerical Modelling of Creep Deformation of Viscoelastic Thick-Walled Cylinder with Fractional Maxwell Model

The deformation of a thick-walled cylinder under pressure is a classic elastic mechanics problem with various engineering applications. In this study, the displacement of a viscoelastic thick-walled cylinder under internal pressure is investigated via analytical as well as numerical modelling. The fractional Maxwell model is initially introduced to describe the creep deformation of high-strength Q460 steel. Subsequently, an analytical solution to the creep deformation of the thick-walled cylinder under both internal and external pressures is deduced with the corresponding principle. The analytical solution is examined with a numerical simulation that incorporates the fractional Maxwell model by a user-defined subroutine. The numerical simulation agrees well with the analytical solution. The limitations of the current study are also discussed.


Introduction
Thick-walled cylinders are widely used in the petrochemical industry, in natural gas, in high-pressure hydraulic systems and in other structures, such as high-pressure oil pipes, petrochemical pressure vessels, heat exchange tubes, storage vessels, nuclear reactor pressure vessels, cannon barrel, steam pipelines and functionally graded materials [1,2]. These internally pressured, thick-walled cylindrical vessels usually operate under hightemperature and high-pressure steam conditions; thus, creep deformation is considered one of the main failure mechanisms of these structures.
The creep deformation of thick-walled cylinders has been intensively investigated. Shinozuka [3] studied the stresses in an incompressible viscoelastic-plastic thick-walled cylinder. Attia et al. [4] studied the residual stresses produced in cast iron cylinders by the creep relaxation of thermal stresses. Schwiebert [5] proposed equations of stress rates in thick-walled cylinders subjected to transient thermal and mechanical loading. Pai [6] used a piecewise linear model to obtain solutions to the steady-state creep of a thick-walled orthotropic cylinder subjected to internal pressure. Sim and Penny [7] used reference stress techniques to investigate the plane strain creep behavior of thick-walled cylinders. Bhatnagar et al. [8][9][10][11] investigated the creep behavior of a thick-walled cylinder under various conditions. Simonian [12] calculated the thermal stresses in thick-walled cylinders, taking into account non-linear creep. Brust and Leis [13] proposed a model for predicting primary creep damage in axial cracked cylinders. Combescure [14] proposed a simplified method for the prediction of creep buckling in cylinders under external pressure. Orlando and Gonçalves [15] derived the lower bound to the creep rupture time of internally pressurized thick cylinders. Most of the above research studies considered thick-walled cylindrical vessels or cylinders with classical creep models. In contrast, few publications can be found which deal with creep in thick-walled cylindrical vessels or cylinders with more powerful fractional modelling techniques.
Thus, in this study, a fractional rheological model is first introduced to describe the creep deformation of Q460 steel under different uniaxial tensile stresses. With this fractional creep model, the creep deformation of a thick-walled cylinder subjected to internal pressure is analyzed through both analytical and numerical modelling.

Materials
Until the first application of fractional calculus in the solution of an integral equation by Abel, it has been a long time since the discussion of fractional calculus between L'Hôpital and Leibniz in 1695 [16,17]. Nowadays, fractional calculus is applied in many science and engineering territories, such as fluid flow, rheology, diffusion, electrical networks, electromagnetic theory and probability. The application of fractional calculus in rheology has been verified to be successful and various fractional rheological models have been proposed to describe the time-dependent phenomena [18]. Among these fractional models, due to its simplicity, the fractional Maxwell model has often been adopted to depict the rheological dynamics for polymers [19], foods [20], fluids [21] and geomaterials [22][23][24]. As fractional modelling is still unfamiliar to non-mathematicians, the constitutive equation for the fractional Maxwell model is presented briefly in the following section.

Constitutive Equations
The traditional integer Maxwell model is composed of a spring element and a Newtonian dashpot in series connection, as illustrated in Figure 1a. The total strain of the Maxwell model includes the instantaneous elastic strain of the Hooke spring and the viscous strain of the Newtonian dashpot. Different from the traditional integer Maxwell model, the fractional Maxwell model [25][26][27][28][29], as shown in Figure 1b, replaces the Newtonian dashpot in the integer model with the Abel dashpot.
in which ε, ε e and ε v are the total strain, elastic strain and viscous strain, respectively; σ, σ e and σ v are the stresses corresponding to the total strain, elastic strain and viscous strain, respectively; α is the fractional order of the Abel dashpot; E is the elastic modulus of the spring; t is the time; η α is the viscosity of the dashpot and the superscript α is utilized to fulfil the principle of dimensional homogeneity.
After simplification, Equation (1) can be rewritten as Equation (2) can be rewritten with the Laplace transform technique as where s, σ(s) and ε(s) are the complex variable, stress and strain in the Laplace domain, respectively. The creep modulus J(s) in the Laplace domain is expressed with stress and strain as follows [25]: Thus, Equations (3) and (4) give Applying the inverse Laplace transform to Equation (5) yields where J(t) is the creep modulus of FMM in the time domain.
If the fractional order in Equation (6) is 1, then Equation (6) reduces to The expression of J(t) in Equation (7) is the creep modulus of the IMM. In other words, the IMM is a special case of the FMM when α equals 1.

Fitting with Fractional Maxwell Model
To verify the practicability of the fractional Maxwell model introduced above, a group of strain-time curves obtained from creep tests of high-strength Q460 steel at various stress levels [30] are digitalized and fitted with the fractional Maxwell model. The creep strain is obtained by multiplying the creep modulus by the constant stress, as shown in Equation (8).
To show the superiority of the fractional Maxwell model, the traditional integer Maxwell model is also utilized to fit the creep curves of Q460 steel, as shown in Figure 2. The fitting parameters are also listed below in Table 1.
From Figure 2, the fractional Maxwell model, illustrated by the blue curve, is found to agree well with the experimental data, illustrated by the red dots, while the traditional Maxwell model, illustrated by the green lines, does not show close consistency with the experimental data. It can also be noted that both the traditional and fractional Maxwell model predict an ever-increasing creep strain; however, there is also a slight difference in that the traditional Maxwell model indicates a constant creep strain rate while the fractional Maxwell model indicates a decreasing creep strain rate. Compared with the traditional Maxwell model, the fractional Maxwell model achieves a better fitting result with one more parameter α; therefore, the fractional Maxwell model seems an appropriate candidate to simulate the viscoelastic deformation of steel.   Figure 3 illustrates a pressured thick-walled cylinder under both internal and external pressures. Under the plane stress condition, the elastic solution of the radial displacement of the cylinder can be expressed as [31]

Methods
where a and b are the inner and outer radius of the cylinder, respectively; p 1 and p 2 are the internal and external pressures, respectively; κ is a parameter depending on plane stress or plane strain conditions; G is the shear modulus and r is the distance from any point within the cylinder to the center of the cylinder. When it comes to the plane strain condition, the radial displacement can be obtained by replacing κ with 3 − 4v in Equation (9) as Hooke's law is decomposed into isotropic and deviatoric parts as where σ m , m are mean normal stress and strain, respectively; s, e are deviatoric stress/strain tensor, respectively; K, G are volume and shear modulus, respectively. The viscoelastic generalization of Hooke's law can be expressed as where the operator D is interpreted as representing partial differentiation with respect to time; f , g, f 1 , g 1 are polynomials in D. The Laplace transforms of the viscoelastic generalization of Hooke's law are Suppose that the solution to a certain problem in elasticity is available; then, the associated viscoelastic solution in which the same stresses are applied at t = 0 to an initially undisturbed body can be obtained by replacing K with g 1 (s)/ f 1 (s) and G with g(s)/ f (s). The actual stresses and strains in the time domain can then be acquired through inverse Laplace transforms. Now, suppose that the material is elastic in hydrostatic compression and follows the fractional Maxwell model in shear. Equation (3) can be rearranged as Comparing Equations (13) and (14), the shear modulus is expressed as Replacing G with Equation (15), the displacement in the Laplacian domain is expressed as wherep 1 andp 2 are Laplace transforms of constant stress p 1 and p 2 . Thus,p 1 andp 2 are expressed as With Equation (17), Equation (16) can be rewritten as Then, the radial displacement in the time domain can be obtained with the inverse Laplace transform of Equation (18) as It can be seen from Equation (19) that the radial displacement is composed of two parts-the elastic part and the time-dependent part. When t → 0, the radial displacement reduces to the elastic component as Equation (20) is almost the same as Equation (10), with only the slight difference that the shear modulus G in Equation (10) is replaced by the elastic modulus E of the fractional Maxwell model in Equation (20), which can be explained by the basic assumption that the time-dependent deformation relies only on the shear stress, as stated previously.

Numerical Implementation of Fractional Maxwell Model
Though the history of fractional calculus is long, it is still somewhat new to nonmathematicians and ordinary engineers. Fortunately, as a general finite element platform, ABAQUS offers users great flexibility to build their own models with a secondary development function. More specifically, when ABAQUS fails to provide an appropriate model to deal with a certain time-dependent material, user subroutine CREEP can be used to define the rheological behavior of this material [32].
If user subroutine CREEP is used to define a material's behavior, the subroutine aims to offer the "uniaxial" creep laws, which will be incorporated into a general time-dependent material formulation. It can be coupled with temperature and electricity analysis. It also should be noted that the meaning of variables defined in this subroutine depends on the specific model chosen. Taking the variables DECRA(1), for example, it may represent the equivalent uniaxial deviatoric/cohesion/compressive creep strain increment when the subroutine CREEP coupled with the metal/capped Drucker-Prager/gasket model. Details of the numerical implementation of the fractional Maxwell model can also be found in a previous study [23] and the manual of the ABAQUS subroutine [32].
To check the reliability of the built-in creep subroutine of the fractional Maxwell model in ABAQUS, numerical creep experiments of high-strength Q460 steel at various stress levels are conducted and compared with the laboratory experiments carried out by Wang [30] et al.
The geometrical model of the steel specimen is illustrated in Figure 4a. The steel cylinder with a diameter 10 mm is 100 mm in length. The steel bar is divided into 100 segments along the axial direction and 20 segments along the circumferential direction, as shown in Figure 4a. The C3D10 is adopted. There are 41,241 elements in total. The steel bar is fixed at the bottom and a tensile stress is applied at the top of the bar, as shown in Figure 4b. The whole simulation is divided into three steps, i.e., initial, elastic and viscoelastic, and the loading process is illustrated as shown in Table 2. In the elastic step, the axial tensile stress increases linearly to the design load and the load will hold constant during the viscoelastic step. Under the stress of 457 MPa, the axial displacement of the steel bar at the end of viscoelastic step is plotted as shown in Figure 4c.  Repeating the simulation process mentioned above, with the mechanical parameters listed in Table 1, the creep of steel bars under different axial stresses is simulated as shown in Figure 5. Form Figure 5, the numerical simulations reproduce the creep experiments of steel bars under different stress levels. Thus, the fractional Maxwell model seems to be a suitable candidate to simulate the first and second stage of the creep deformation of steel bars.

Numerical Simulation of Creep Deformation of Thick-Walled Cylinder under Internal Pressure
A plain strain model is established to simulate the creep deformation of a thick-walled cylinder under internal pressure. The inner and outer radii of the cylinder are 1 cm and 6 cm, respectively. The elastic modulus of cylinder is 10 GPa and the Poisson's ratio is 0.3. The viscosity of the cylinder is 10 13 Pa · s a and the fractional order is 0.5. Each edge of the cylinder is divided into 50 segments. A bias mesh scheme is adopted on two straight edges and the bias ratio is 10, as shown in Figure 6a. The internal and external pressures with the magnitude 10 MPa and 5 MPa are applied on the inner and outer surface, respectively, and a symmetric condition is cast on the two straight edges, as shown in Figure 6b. The internal pressure is held for 3600 s. During the whole loading process, the radial displacement at the inner surface of the cylinder is recorded. The radial displacement at the end of loading is plotted as in Figure 6c. It should be noted that Figure 6c is plotted in a local cylindrical coordinate system. The theoretical radial displacement at the inner surface of the cylinder is plotted according to Equation (19) and radial displacement by the numerical simulation is also plotted with the data obtained at each step, and they are both plotted and compared in Figure 7. In Figure 7, the numerical simulation results almost coincide with the analytical solution predicted by Equation (19). The analytical solution of the displacement at the inner surface under constant pressure is validated by the numerical simulation.

Discussion
In line with the previous section, it should be noted that there are also limitations of the current study. Firstly, only the viscoelastic deformation is considered in this study and the plastic deformation is not taken into consideration. This is a requirement of the elastic-viscoelastic correspondence principle [33]. This is also a requirement of common practice for most containers and pipes that operate below the yield stress. Secondly, only the transient and the steady-state creep stages are considered both in the analytical solution and numerical simulation, and the acceleration creep stage is omitted for simplicity. Both the traditional integer Maxwell model and the fractional Maxwell model can only describe the transient and the steady-state creep stages; more complicated rheological models may be introduced to tackle the complexity involved with the acceleration creep stage. The numerical simulation with more complicated rheological models may not be overly difficult to implement; however, whether there would be analytical solutions generated by these models remains uncertain. The fractional Maxwell model seems to find a balance between more complicated models and practical needs, as the fractional Maxwell model captures the primary features of the steel creep and also enables an analytical solution of the displacement for the internally pressured, thick-walled cylinder.

Conclusions
The viscoelastic deformation of a thick-walled cylinder under internal pressure is investigated via analytical modelling and numerical simulation. Based on the study, the following conclusions are made:  Data Availability Statement: Data sharing is not applicable to this article.